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Time Dependent Scattering from a Grating 
Li fan 1 , Peter Monk 1 


Abstract 

Computing the electromagnetic field due to a periodic grating is critical for assessing the perfor¬ 
mance of thin film solar voltaic devices. In this paper we investigate the computation of these fields 
in the time domain (similar problems also arise in simulating antennas). Assuming a translation 
invariant periodic grating this reduces to solving the wave equation in a periodic domain. Mate¬ 
rials used in practical devices have frequency dependent coefficients, and we provide a first proof 
of existence and uniqueness for a general class of such materials. Using Convolution Quadrature 
we can then prove time stepping error estimates. We end with some preliminary numerical results 
that demonstrate the convergence and stability of the scheme. 
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1. Introduction 

We wish to approximate time domain electromagnetic scattering from a periodic grating. We 
shall assume that the grating is translation invariant in one direction, so that Maxwell’s equations 
can simplified to obtain two Helmholtz equations governing the s- and p-polarized waves. Our 
intended application is to modeling solar voltaic devices. Usually such devices are modeled in the 
frequency domain, and for descriptions of frequency domain applications in this area see [1, 2]. We 
will consider the problem in the time domain with the potential benefit of being able to compute 
results at a range of frequencies in one simulation although we do not investigate that aspect 
here. Note that although our interest is in periodic gratings, similar problems also arise in antenna 
theory [3] and [2, SectionlO.2.2]. We expect that the theory developed here can be extended to 
that case. 
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We start by describing the problem assuming that the grating is translation invariant parallel to 
the y axis, so that the permittivity of the material in the grating is independent of y (the magnetic 
permeability is assumed to be that of free space). Because of the reduction in dimension afforded 
by translation invariance, Maxwell’s equations can be reduced to the wave equation in two spatial 
dimensions (see for example [4, Section 5.1]). So we assume that the electromagnetic wave is 
described by a scalar function u = t/(x, t) depending on position x = (x, z) E M 2 and time t > 0 
that satisfies 


1 O^u 

^2 b *^ = V -( a *Vu)m 


t > 0. 


( 1 ) 


c* dt 2 

Here c > 0 is the speed of light in vacuum, the symbol * denotes convolution in time and the 
functions a and b describe the medium in which the electromagnetic field propagates. Two choices 
are of interest: 


1. The choice 

a = 5(t) and b — e r 

where 5 is the Dirac delta and e r is the time domain relative permittivity of the medium. In 
this case the wave is said to be s-polarized and u represents the y component of the electric 
field. 

2. Alternatively 

a — l/e r and b — S(t ) 

in which case the field is said to p-polarized and u represent the y component of the magnetic 
field. 


Often, for simplicity, it is assumed that a and b are independent of frequency and are real, bounded 
and uniformly positive piecewise continuously differentiable functions of position. However for 
realistic materials both a = a(x,t) and b = 6(x, t). 

Since the medium is assumed to be a grating, there is a period L > 0 such that 


a(x + L, z, t) 
b(x + L, z, t) 


a(x, z,t) 1 x — ( x z \ ^ ]^ 2 an q £ E M. 

b{x,z,t) J 


In addition the grating is assumed to have a finite height H such that 


a — b — 5 for all x E R and z < 0 or z > H. (2) 

This assumption can easily be relaxed to allow for different materials above and below the cell (one 
of our examples features this). We postpone further discussion of the assumptions regarding these 
coefficients until we have introduced sufficient notation. 


Later we will reduce the problem to a bounded domain called the “unit cell” defined by 

Q = (0, L) x (0,if), 

and we will also need the unbounded strip 

S = (0, L) xM. 

We assume that the total field u is due to an incident plane wave u l propagating towards the 
bottom y — 0 of the grating. In particular 

7i z (x, t) = f{t — d • x/c), x E M 2 , 


2 


for some twice continuously differentiable function /, and unit vector 

d = (di, cfo) :=(cosscq sin<a), 0 < a < tt. 

So = 7t/ 2 gives a plane wave propagating up along the z-axis at normal incidence to the grating. 
In general, by linearity, the total field can be written as 

u = u l + u s 

where u s is an unknown scattered field to be determined. 

Notice that the incident field u l is not periodic in x but instead 

u l (pc + L, z, t ) = f(t — di(x + L)/c — d 2 z/c) = f(t — d\L/c — d • x/c) 

= u l (x, z,t — diL/c) (3) 

for any x, z and t. Thus we expect that the scattered field and hence the total field u should have 
the same translation properties so we impose 

u{pc + L, z, t) — u(x, z,t — d\Ljc ) for all x, z and t. 

This is the time domain counterpart of quasi-periodicity in the frequency domain [5] . 

To avoid the retarded periodicity condition (3) it is common [6, 7] to change variables for this 
problem and define 

w(x, z, t) = u(x, z,t + (x — V)d\jc ), and w l (x, z, t) = u l (x , z,t + (x — L)d\/c). 

With this definition it is clear that w and w l (and hence w s defined in the same way) is periodic 
in x since 


w(x + L, z, t) = u(x + L, z, t + xdi/c) = u(x, z,t + {x — L)d\/c ) 
— w(x,z,t). 


The incident field becomes the trivialy x-periodic field (independent of x) 

w l (x, z, t) = u l (x, z,t + (x — L)di/c) = f(t — Ldi/c — d^zjc). (4) 

We assume the field in the unit cell is quiescent before t — 0 so that w s — 0 for t < 0. This requires 
that w l — 0 on ft for t < 0 or 

f(t — Ldi/c — d2z/c) = 0 

for t < 0 and 0 < z < H. Assuming d\ > 0 and g ?2 > 0 it suffices that f(t) = 0 for t < 0. 

With the above change of variables equation (1) become more complicated. Setting ei = (1, 0) we 
have 


_ _ d\ dw 

vu = Vw ---ei, 

c dt 


A u = 


. d\ d 2 w d? d 2 w 

Aw-2 — ^^ + 1 


c dtdx ' c 2 dt 2 
So since V • a\7u = aAu + Va • we conclude that w satisfies 

b — ad?\ d 2 w d\ d ( dw\ 

1 * = V • (a * vw) -— a * 

c ox v 1 


dt 2 


dt J 


d\ d 2 w 
~~ ^ a *~dtdx 


( 5 ) 
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for all x, z and t. The problem we wish to solve is to compute w such that (5) holds together with 
the splitting 

w = w l + w s in S x M, 

together with the initial conditions 


dw s 

w s = —— = 0 in S at t = 0, 
at 

and in addition w is L-periodic in x. 

The change of variables we have used above is well-known in the engineering literature [6] where an 
explicit finite difference time domain (FDTD) technique is used for discretization and time stepping. 
In this case the time step must be chosen depending on the angle of incidence. Some later authors 
[3, 8] have use implicit schemes to avoid stability restrictions. We shall use implicit methods here. 
Another interesting FDTD paper for gratings that discusses, amongst other things, locally refined 
grids is [9]. The numerical analysis of several models of dispersive media are discussed in [10, 11] 
and in particular [12], but these studies do not cover grating problems, and do not provide an 
analysis for a general class of problems. 

Existence and uniqueness questions for gratings were studied in [ ]. In that thesis, existence and 
uniqueness is proved using a weighted norm in time for frequency independent coefficients. This is a 
special case of our result which improves the norm and also allows frequency dependent coefficients. 
We use the Laplace transform as a tool and in addition prove error estimates for a limited class 
of time-stepping schemes. We also provide preliminary numerical results verifying the temporal 
convergence rate of the scheme, as well as showing results for two frequency dependent models 
covered by our theory. 

The layout of the paper is as follows. In the next section we prove existence and uniqueness 
of a solution to the time domain problem. We also give error estimates for Backward Euler and 
Backward Differentiation Formula 2 (BDF2) based discretization in time. Both these time stepping 
rules are implicit, the former being first order and the latter second order [13]. In Section 3 we 
show how we discretized in space and truncated the problem using a domain decomposition strategy. 
In Section 4 we give three numerical examples: we first verify the predicted convergence rate in 
a simple case, then we give two examples of computations using standard frequency dependent 
coefficient models. Finally in Section 5 we put the study into perspective. 

To simplify the presentation, for the remainder of this paper we will assume s-polarization so that 
a = 5. The p-polarization case can be analyzed in a similar way. In that case b = 5 and we need to 
required that 1 — is strictly positive. In addition, we need to ensure coercivity of the appropriate 
bilinear form. The main assumption is that if a(s)is the Laplace transform of a with parameter s 
(see upcoming (11)) and if a = 5R(s) > 0 then there is a constant C depending on a such that 

»(sa(s)) > C. 

The corresponding assumption for s-polarization is given in detail Assumption 1, and we would 
also need to assume the differentiability conditions on a. Numerical examples of the p-polarized 
case are not investigated here. 
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2. Existence and uniqueness 


For theoretical purposes it is convenient to work with the scattered field. So setting 

w — w l + w s 

and noting that w l satisfies (5) with a = b = 5 we have that w s satisfies 


where 


b — d 2 5^ d 2 w s 
dt 2 
w s 
dw s 
~dt 

w s (L , z, t ) 
dw s 


F = 


= A w 


d\ d 2 w s 


c dtdx 
— 0 in S at t — 0 

= 0 in S at t = 0 


+ F in S x 


w s ( 0, z, t ) for t > 0, z G R, 
Qw s 

— —(0,2:, t) for t > 0, z G R, 
ox 


5(t ) — d 2 w l 

* 


( 6 ) 

( 7 ) 

( 8 ) 
( 9 ) 

( 10 ) 


c* J dt 2 

Notice that, by our definition of H (see (2)), F = 0 if z > H or 2:<0 so F has support in ft. Our 
assumptions on / also imply that provided b is causal, F — 0 in ft for t < 0. We can extend w s by 
zero to time t < 0 and hence obtain a causal function defined for all t. 

To analyze this problem we use the Laplace transform in time [14, 15]. For any sufficiently smooth 
function g = g(t) with at most exponential growth for large time, the Laplace transform is g — g(s ) 
given by 

r oo 

g(s)=C(g)(s)= g(t) exp(-st) dt, (11) 

Jo 

where we choose the transform variable to be s = o — iuo for o G R and o > 0, and uo G R. Then 
because of our initial conditions 

dw s 


' d 2 w s \ 2 / d 2 w s \ 

c[ ^fi)= ,w a " d£ («ai) 


dx 


In addition 


n oo 

w l (x,z) — / exp(— st)f(t — Ldi/c — d 2 z/c) ds 

Jo 

/ oo 

exp(— s(r + Ldi/c + d%z/c)) f (r) dr 

-Ldi/c—d2z/c 

= f(s) exp(—sLdi/c) exp(—sd 2 z/c). 

As expected this is a scaled plane wave in the Laplace domain that is actually independent of x 
(see equation (4)). 

Formally taking the Laplace transform of (6) we are lead to seek w s G H^(S) that satisfies 


-dl\ , s 

cf) W 

= Aw® 2**“’ +FinS, 
c dx 

w s (L,z) 

= ii) s (0, z) for y G R, 

dw s 

dx ™ 

dw s 

= ~—(0 ,z) for £ € R. 

dx 
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Here F = s 2 (1 — b)w l / c 2 and b is the Laplace transform of b. 
To formulate a variational problem for this system, let 


H^S) = {feH 1 (S) | f(L,z) = f(0,z) for all zeR}. 

(where the subscript p recalls the x-periodicity of the functions) with the 5 -dependent norm 


\w 


I H$(S) 


|Vu ) s | 2 + 



dA 


where dA — dxdz. We shall require that w s E H^{S) and this requirement replaces a radiation 
condition. It holds because = a > 0. 

We can now write a Galerkin formulation for the Laplace transformed problem as usual by multi¬ 
plying by the complex conjugate of a test function and integrating by parts (using the periodicity 
of the normal derivative to cancel terms on x = 0 and x — L). We seek w s G H^(S) such that 


a(w s ,C = f F£dA for all f € HUS), 

Js 

where the over-bar denotes complex conjugation and 



Vw s • V£ + s 2 



w s f + 2 s 


d\ dw s 
c dx 


€ 


dA. 


( 12 ) 


At this stage we need to specify our assumptions on the frequency dependent coefficient b. We 
assume 

Assumption 1 . The coefficient b = 6 (x, s) is piecewise continuously differentiable in x. In addition 

1. For almost every x G S the coefficient b is analytic in s for $ft(s) > cro > 0 for any ao, and 
bounded independent of s (but the bound may depend on cfq). 

2. There is a constant 70 such that 

3?(s(6(x, s) — d 2 )) > cryo > 0 

for $ft(s) = a > 0 and all x G S. 

3. b(x, z,s) = 1 for z < 0 or z > H and all x, s. 

In Section 4.2 we will give two important examples of a frequency dependent coefficient satisfying 
the above assumptions. 

Now we can use the Lax-Milgram Lemma to prove existence of a solution. First we verify coercivity 
and continuity. 

Lemma 1 . Suppose that b satisfies the Assumption 1, and 5ft(s) = a > > 0 for some ctq. Then 

the sesquilinear form a(-,s-) is coercive and bounded. In particular for every v G H^(S) 

|a(D,sD)| > crmin(l,7o)||u||^i (s) . 

In addition there is a constant C depending on <7o but independent of s, u,v E H^{S) such that 

|a(u,D)| < C , ||u|| if i (5) ||u|| H i (s) . 
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Remark 1 . The dependence of the constant in the boundedness estimate above on ctq comes from 
the frequency dependence of b. If b is frequency independent, then the constant is independent of 
cr 0 . 


Proof. Using the by now standard trick of Bamberger and Ha Duong [ 5] we choose £ = sv and 
obtain 


a(v , sv) = 


L 


s |Vh| 2 + s|s| 


b ~ d l \ I~|2 , oi 1 2 d l 9v — 

-o M + 2 s — ^-V 

c z c ox 


dA. 


But since integration by parts in x shows that 


L 


dv- „ dv , . 

—v + v— dA = 0 

g OX ox 


we have 



and so we have proved coercivity because, using our assumption on 6, 


5ft[a(i), sv)] = [ 

Js 


5ft 


s |Vf)| 2 + \s\ 2 s 


b-di 


dA 


> o min(l, 7 o)||' 0 ||^i (s) . 


In addition a(-, •) is bounded because for any u,v G H^(S) we have 

Isl 2 * 

|a(u,t))| < ||Vu|| L 2 (s) ||Vu|| £ 2 (s) + -fr\\ b ~ d?|U~(S)||h|| L 2 (s) ||'i)|| i 2 (s) 

+2|s|-i||Vu|| L 2(5)||i)|| i 2 ( S ). 


Use of our assumption on the coefficient b and standard inequalities demonstrates the required 
continuity. □ 

Using the above lemma and the Lax-Milgram Lemma we have the following result 

Theorem 1 . Assume that b satisfies Assumption 1 and 5ft(s) = o > 0. Then problem (12) has a 

unique solution and 

amin(l,7o)|K|| H i ( 5' ) < ^ II& “ ilU-iS) Ill’llL 2 (5) 


This verifies the existence of the solution in the Laplace domain and this in turn gives a time domain 
existence and uniqueness theorem. Note that we need the condition 5ft(s(6(x) — d\)) > cryo > 0 for 
all x in S. Since d\ = sin# where 6 is the angle of incidence, this may limit the magnitude of the 
angle of incidence. 

Using Lubich’s theory of convolution quadrature [14], and considering a finite time period [0,T] for 
some T > 0 we can now state the following existence and uniqueness result. The usual space for 
this theory is 


H^((0,T),Hl(S)) = {g| (0 ,T) I g e with g = 0 for t < 0}, 

which places compatibility conditions on the data at t = 0 if m > 1. 
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Theorem 2. Suppose F G iL™((0, T), L 2 (*S)) for some m > 0. Then there exists a unique weak 
solution w s G ^^((O, T), H^(S)) of the time domain problem (6)-(10). 

Proof. This is an application of a slight generalization of [14, Lemma 2.1] (for similar analysis see 
Section 2.2 of the same paper). In the transform domain we can write 

w s = K(s)F 

and the solution operator K(s) : L 2 (Q) H^(Q) satisfies, for any fixed cfq > 0 and all s with 

9?(s) > cto 

ll^( 5 )l|L 2 :if 1 < C 

for some C depending on ao but independent of s. Here || • is the operator norm for maps 

from L 2 (S ) -G Hp(S). Parseval’s theorem gives, using the contour = a > cro, 

|| exp(-at)w s \\ H rn^ 0 ^ H i^ S )) < C\\ exp(-at)F\\ H m mT ^ L 2 {s)y 


□ 


To obtain results for time discretization we can appeal to Lubich’s theory of Convolution Quadra¬ 
ture [14]. This is simplest for a multistep scheme which we now describe (an implicit Runge-Kutta 
scheme could also be used but is beyond the scope of the paper). Let At > 0 denote the timestep 
and let t n = nAt, n > 0. Given g(t , y), consider the differential equation y' — y(t, y ) for t > 0 with 
y(0) = 0, then a general fc-step multistep scheme applied to this ordinary differential equation is to 
find {y n }^Lo such that 

k k 

a jVn-j = At / 3jg(t n -j , y n -j) 

3=0 3=0 

where we assume compatible initial data and set yj = 0 for j < 0. The coefficients {cvy, /5y}y =0 
define the method and we assume that <ao/A) > 0. Then following Lubich [14], let 


7(C) = 


Y.UW 


CgC. 


Lubich shows Convolution Quadrature time stepping is equivalent to the following parameterized 
problem. Let 

W s (x)=^<(x)C 

3=0 

for CgC small enough. Here we shall prove that converges to w s (-,t n ) as At 0, and we take 

— 0 for n < 0. Then W s G H^(S) satisfies 


a A t(W s ,0 = J dA for all £ € flJ(S), 


(13) 


for all |Q < 1, where 

a A t(W s ,£) = [ 

Js 


vr • v£ + 


( 7(0 y ( bAt - d 

V At 


W s £ 


7(0 <h dW* ? 
At c dx ? 


gL4, 








and b/\t and F^t are obtained from b and F by replacing the transform parameter s by 7 (C)/ At. 

To obtain a time stepping scheme it then suffices to equate terms in (/ n on the left and right 
hand side of (13) to obtain equations for w * in terms of previous values. When b is frequency 
independent, we obtain the usual multistep update formula that can be derived alternatively by 
applying the multistep method to the standard weak form of ( 6 ) (however the approach we have 
given will provide an error analysis). 

Obviously this scheme is not yet computable because S is unbounded. One possibility is to rewrite 
(13) by truncating S using auxiliary boundaries above and below the grating, then using an integral 
equation or Dirichlet-to-Neumann map on the auxiliary boundary to close the system. We discuss 
an approach similar to the Dirichet-to-Neumann map approach based on upwind transmission 
conditions in the next section. 

The theory of convolution quadrature provides an error estimate for this problem. Let (d/\tW s ) n 
denote the coefficient of ( n in the expansion of ^(QW S /At. This is the finite difference approxi¬ 
mation of the time derivative corresponding to the given multistep method. We have the following 
theorem: 

Theorem 3. Suppose the multi-step method is A-stable, order p convergent and that 7 (C) has no 
poles on the unit circle. In addition suppose b satisfies Assumption 1, that m = p + 2 and that 
d^F/dP =0 at t — 0 for j — 0,1, • • • , m — 1. Then for 0 < t n < T we have 

11 v(» s (,g - oilmen) + l|i|» s (.,g - ( d At w s ) n )\\ LHQ) 

[t Qmp 

< c(a ty> f \\— (.,T)|| L 2 (n) dT 

where C is independent of At, w s and the discrete solution, but depends on T. 

Remark 2. The canonical examples of suitable time stepping schemes that satisfy the requirements 
of the theorem are backward Euler (p — 1) and BDF2 (p — 2). In the latter case m — 4 and we 
need F to have 4 weak derivatives in time and satisfy the compatibility conditions F/dP — 0 at 
t — 0 for j = 0,1, ■ ■ ■ ,3. 

Proof. The proof is a direct consequence of Theorem 3.1 (see also the discussion before Corollary 
4.2) of [14]. □ 

3. Discretization in space 

We can derive a time stepping scheme by truncating (13) using integral equations or the Dirichlet- 
to-Neumann (DtN) map and then equating terms in £ n . However in order to demonstrate the 
viability of the approach we will use the Laplace domain approach of Banjai and Sauter [16]. 
This requires to solve the Laplace domain problem for several different choices of the transform 
parameter s. The algorithm is exactly as in Banjai and Sauter’s paper and so we only give details 
of the Laplace domain problem that we solve. As pointed out by Banjai and Sauter, this approach 
is trivially parallelizable, although we have not investigated this or other algorithmic improvements 
here. This can help improve the elapsed time for running the algorithm. 

First we make a minor change and will solve for the total field in D (this avoids evaluating F 
everywhere in the grating, but is less convenient for analysis). Let c denote the speed of light below 
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the grating for z < 0. The x periodic total field w = w(x, z) satisfies 


Aw — 2-d\w x — %(b — d\)w = 0 in 11, 
c c z 


(14) 


To accommodate one of our examples, we assume a slight generalization of the problem discussed 
so far. We assume that b = 1 for z < 0 and b = b\ for z > H where b\ is spatially constant but 
could be frequency dependent and so depend on s. In the region 11 we have, in general, 6^1. 

For z<0 we see that w satisfies 


Aw — 2-d\w x -^ (1 — d\)w — 0. 


(15) 


This field can be decomposed into an incident field given by 

w l (x , z) = exp(— s—z) exp(— s—^) f(s ) 
c c 

and scattered field as before so that w = w l + w s . The scattered field also solves the above 
differential equation for z < 0 or z > H and we now derive an equation for the scattered field for 
z < 0. Since w s is periodic 

w s {x,z) = tin( z ) exp(z27mx/L) 
ne N 

for suitable Fourier coefficients Substituting this expansion into (15) gives the condi¬ 

tion 

d? 


dz* W * ° 


- + ^2n7r/L + ^ ti s n = 0, y < 0. 


Let K s n be defined by 


«) 2 = (^) (l + (2rwrc/ (Ls) + idi) 2 ) 


We need to choose the signs of the square root so that we have a decaying solution as z —» — oo. 
Then 

w s (x,z) = y]<exp(i27mx/L)exp«z), for z <0. (16) 

nG N 

for suitable expansion coefficients w!^ and where $ft(/^) > 0. 

For z> H, the transmitted total wave w f satisfies 


A w t — 2-d\w i L — %(bi — dVjw 1 — 0 
c c z 


(17) 


There is no incident wave and the wave is only transmitted. Proceeding as before we define W, 
by 


( 4> 2 = © lfc + 1 


( 2n7rc 


Ls 


+ id\ 


where we choose again 9?(/^) > 0. The transmitted field for z > H is given by 

w f (x, z) = exp(i27rnx/L) exp (—^ n {z — H)), for z > H > 0 

ne N 


(18) 
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for suitable expansion coefficients w^. 

The fields w s for 2 ; < 0, w 1 for z > H and w for 0 < z < H are related by continuity conditions 
across the interfaces at 2 : = 0 and y — H. At z = 0 we need 

w(x, 0) = w l {x , 0) + w s (x, 0) and w z (x, 0) = w l z (x, 0) + w s z (x, 0), 0 < x < L. 

At the upper interface z — H we need 

w(x, H ) = w f (x, H ) and w z (x, H ) = w l z [x^ iL), 0 < x < L. 


To formulate a problem suitable for spatial discretization we have adopted the one-way wave equa¬ 
tion or characteristic equation approach of [17] which is based on impedance boundary conditions 
(see also the UWVF [18] and for strongly related methods [ 9, 20]). 

Let y be a positive constant at our disposal (in fact we choose y = 1). Given Aq G L 2 (0, L) and 
G L 2 (0, L) define w = w( Aq , \^) G iLp(fl) by 

5 5 2 

Aw — 2-d\w x -2 (b — d\)w = 0 in Q. 

dw s A _ _ TT 

— - 1 — yw = \tt lor z = H. 0 < x < L. 

dz c 

dw s 

— — -b -yw = A7 for z = 0. 0 < x < L. 

dz c 

This is easily written as the variational problem of seeking w G iLp(fl) such that 


B(w,v) — f(v) for all v G iLp(fl) (19) 

where 



Here £# = {{x, H ) | 0 7 x 7 -l) while — ((x^ 0) | 0 7 x 7 1/) , and 



X H v ds + 



A 0 v ds. 


For z > H define ^(A^) by requiring that 

dw 1 Q 

— - rjw 1 — At, for z = H, 0 < x < L, 

dz c 

together with the expansion (18). Using a trigonometric expansion (noting the periodicity of the 
solution and its normal derivative) 


A b = E A 7 n exp(i 27 rnx/L) 
ne N 
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then the Fourier coefficients of the transmitted field are 


U>n = ~ X H, n /( K n + s V/c), HE Z, 


and we can see that this field is aways well defined since lifts > 0 and 5ft > 0. Define Fh{ A^) 
by 


Fh{\ + h ) 


dw t s w 

-1- -r\w on E# 

oz c 


T ^ H,n 

nG N 


( K n - (sfc)y \ 

V K n + ( s / c )vJ 


exp(i27rnx / L ). 


Similarly for z < 0 define w s ( Aq") by requiring that 

dw& g 

—-- rjw s — Art for z = 0. 0 < x < L. 

oz c 

together with the expansion (16). Suppose 

A » =E A+ n exp(i 27 rnx/L) 

nG N 

then the nth Fourier coefficient of the scattered field is 


Define F 0 ( Aj) by 


n(A+) 


K = -Aj n /« + (s/c)ri). 


dw 


_ + -T]W S on H 0 

oz c 


E A t 

nG N 




(s/ 0)7] 


K n + (s/c)7] 


exp(i27rnx/L). 


It remains to derive equations for A^ and AThis is done by enforcing the transmission conditions 
in impedance form. At z = H we require 


dw 

dz 

dw 

dz 


s ^ 

H —rjw 
c 

s A 

- rjw 

c 


dw 1 s 

— -1 —riw at z — H, 

dz c 

dw 1 s , 

— - rjw at z — H. 

dz c 


Writing these equations in terms of the unknown functions, 


g 

A# ~ 2-77u)(A 0 , A^) — A^, 

^h — Fh(^h), 


where we have avoided computing the normal derivative of u by writing 


dw 

dz 
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The same process can be applied at z = 0. We require that 


dw 

——b ikrtu 
oz 

dw ^ 

— - ikrjw 

oz 


d(w s + w l ) 
dz 

d(w s + w l ) 
dz 


+ ikrj(w s + w l ) at y — 0, 
— ikrj(w s + w l ) at y = 0. 


Let 

dw 1 s ^ dw 1 

j- — — -1 —rjw and /+ = — - rjw at z = 0 

oz c oz c 

This gives the equations 


-A 0 - = —F 0 (A+)-/_, 

— A 0 + 2-t]w(\q , \ H ) = — Aq“ — /+. 

In summary we must find A^ G ^2(^77) and A^ G £2(^0) such that 

A# — 2-77ih(A 0 , A^) = A^, 

= 

A 0 " = F 0 (A+) + /_, 
A 0 — 2-rjw(\ 0 , \ H ) = Aq" + /+. 


In order to discretize the problem, we expand all boundary functions as a finite Fourier series: 

N 

X h’ N = ^f,n eX P(* 27r W L )> 

n=—N 

N 

X o’ N = A L ex P(i2*nx/L), 

n=—N 

and for ease of notation define ^{x) — exp(i 27 rnx / L)\^ H and ^{x) = exp(z27rna;/L)|s 0 . Then 
using the inner products 


{u,v)h 



uv ds , 



uv ds , 


we have the discrete system 

{ x ~h N X)h - 2 ^v(™( X 0 ,N ’ X h N )^p)h - ( A h N ^?)h 

(\ h ’ N ,iP?)h-(Fh(^ N ),^)h 

(Ao’ iV ,^ r 0 )o-(Fo(A+’ iv ),^ 0 )o 

(Aq j ^3)0 — Aq , X H ),rpg)o — (Aq", )o 


0, — iV < p < N, 

0, ~N <q<N, 

</-,$} 0 , - N<r<N, 
(/+,€) 0 , ~N<s<N. 


In our calculations we use a finite element approximation to u; based on standard quadrilateral 
elements and continuous mapped piecewise bilinear functions. The above system is solved using 
deaI.II [21]. 


In the upcoming calculations we can choose N to be relatively small (in fact N = 10 in the first 
experiment) because we expect only a few propagating modes in the frequency domain. Then N 
needs to be chosen to include these modes and a few evanescent modes in addition. 
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4. Numerical Results 


4-1. Quantification of the error 


We start with a simple problem with a known exact solution to test convergence. We fix c and 
suppose e = 6-5(t) for z < hi and e = 6+5 (t) for z > h{ for some fixed interface height hi with 
0 < hi < H. Consider an incident field at normal incidence (so d\ — 0) defined by 

u l (x,z,t) = f - hi)+t^j , z < hi, 

where / is a given function. Then the scattered field for z < hi will be 

u s (x,z,t ) = g s — hi) + t 

for some function g s and the transmitted wave in z > hi is 

u\x, z,t ) = g t ( - ~y - hi) +t 

for some function gt. Continuity of the fields at z = hi implies 

fit) + g s (t) = g t (t) for all t. (20) 


Continuity of the normal derivative implies 


—/'(*)-—■ 


■9s(t) = 


Hence, integrating this 







9t(t) 


-\/e+iV) + V^+9s(t) = -y/eZg t (t). 


Using (20) 


9s (t) = 


v^+ + 


fit) 


and 

9t(t) = f(t)+g a (t) = f{t). 

V e + + v e - 

We choose f(t) to be sufficiently smooth for our convergence theorem to hold and to vanish for 
t< 0, and in particular for parameters m, a[ nc > 0 and /3- inc we define 


/(*) = 


0 for t fiinc Or t 7T f C^inc "b /^inc 

sin (c^inc(^ Anc)) for Anc t 7T / 0^j nc H - (3 

inc 


( 21 ) 


This function is in ib m (M), and we choose m — 4. We choose a[ nc = 4 and Anc = 0.5 in this 
example. 


We can now solve the above problem with c = 1, e_ = 1 and 6+ = 4 to generate a solution to the 
scattering problem and hence test convergence of the time stepping method using the fixed spatial 
mesh in Fig. 1 and BDF2 in time. The final time is T = 4 by which time the wave has essentially 
exited the computational region. Results are shown in Fig. 2. Both in the L? and H 1 norm the 
convergence rate is ultimately 0((At) 2 ) as expected from Theorem 3. No instability is evident. 


14 










Figure 1: Simple almost uniform spatial mesh used to investigate convergence of the time stepping 
scheme. The domain is D = [0,1] x [0,1], and the interface is at hi = 1/2. 


4-2. Frequency dependent materials 


Our Assumption 1 on the coefficients in the differential equation handles at least two important 
cases relevant to practical applications to solar-voltaic components. To see why this is necessary 
note that in thin film devices the size of components is close to the wavelength of light. At these 
frequencies (for example corresponding to a free space wavelength of 500nm) metals can no-longer 
be modeled as perfect conductors and the light penetrates an appreciable distance into the metal. 
For example, at a free space wavelength of 500nm, e r — —2.4683 + 3.1173i for gold. So in fact e r 
is often complex valued and the real part may not be positive. In a Drude model (commonly used 
to model metals [ 22 , Ch. 2 ]) we have 


(s) — CKm T 


Pn 


S(1 +7mS) 


( 22 ) 


for positive, perhaps spatially dependent, real constants a m , /3 m and y m . Note that if = 0 this 
reduces to the usual model of conductivity. To verify property 2 of Assumption 1 in the case of a 
Drude model note that 


$ft(s( 6 m (s) - d\)) = 3? ( sa m + 


(1 + 7 m s) 
|1 + 7ms| 2 


— dis 


= a - d 1 + 
> o' (a m - d\) . 


(j| 1 + 'JmS \‘ 


r(l + CT7 m) 


Provided, for example, a m — d\> 70 > 0 we have the desired positivity. 


Because of the wide range of frequencies need to simulate a solar cell components across the solar 
spectrum, it is also necessary to take into account frequency dependence for dielectics. Dielectric 
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Figure 2 : Time step convergence of the method. Top row: The absolute L 2 (left) and H 1 norm 
(right) error at different times. Bottom row: to check the convergence rate we show a log-log plot of 
the absolute L 2 (left) and H 1 error (right) as a function of At at t = 1.5. We also show a 0((At) 2 ) 
reference line. 


components are often modeled as having no absorption but frequency dependence (so $s(b(s)) = 0 ). 
A commonly used model is the Sellmeier equations [23, page 472]. In this case, the simplest model 
is 


b s (s) — 1 + 


cr © 


(23) 


1 + (3 s s 2 ' 

More generally there are usually sums of rational functions of the same form as above. Here a 
and P s are positive real constants. This model also fits into the theory because 


3?(s(S s (s) — d \)) = §? ( s + a 


= a 1 + a 


(s + ^|s| 2 s) 
|l + /3 s s 2 | 2 
(l + /3 s |s| 2 ) 


— sd\ 


-d\ 


|1 + /3 s s 2 | 2 

> cr(l — d\). 

So provided 1 — d\ > 70 > 0 we again have the necessary lower bound. 


Neither of these models satisfy the Kramers-Kronig relationship that guarantees stability and 
causality, but, as we have proved, they still provide a stable time dependent response for finite 
time. Note that the Cauchy model of a dielectric [23, page 468] does not fit into this theory. 
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Figure 3: Two gratings used for the numerical results on frequency dependent media. In this figure 
ni, 7225 na and Tig label the refractive indices for the various subdomains. The permittivity is 
obtained by squaring the refractive index so e r = tt^ in the grating in the left example. The left 
figure shows a simple metallic grating in air. The right figure shows a cylinder slightly above a 
glass substrate and features a frequency dependent refractive index in the infinite glass region. 

4.2.1. Scattering from a Drude metal 

We now show an example of a typical grating (without other components of a thin film solar cell 
in order to emphasize grating effects). The geometry of the experiment is shown in the left panel 
of Fig. 3. Light is incident (at incidence 6 — 6°) on the thin metal grating which is modeled by 
a fictional Drude medium (see (22)) having parameters a m = 4.0, /3 m = 10.0 and = 0.5. We 
set 77/2 and rig to be 1 and chose 10 modes for the top and bottom boundaries. By setting 777/ = 4, 
cq nc = 4.0 and /3- mc = 0, we chose (21) as the incident function for this example. In particular 
the incident field has a non zero Fourier transform at low frequencies where our fictitious Drude 
model has negative real part, so that the dispersive and dissipative nature of the medium is probed. 


The domain is D = [0,1] x [0,1] and so, using the notation in Fig. 3 left panel, L — 1. In addition 
we choose L\ — 0.1, L m = 0.1 and L g — 0.05. The speed of light is set to c = 1 and we integrate 
until T — 4 using 512 time steps so At = 0.0078. The mesh is from Fig. 4 left panel. This spatial 
mesh was generated using gmsh [24] a triangular mesh generator that can post-process the mesh to 
create a quadrilateral mesh. Obviously the resulting mesh is rather poor but this is useful to test 
the sensitivity of the method to mesh perturbations. Density plots of the computed total field are 
shown in Fig. 5. At early times the incident field is clearly visible followed by a strong scattered 
field (in a thin film solar cell, other structures would serve to trap the energy). In addition the 
transmitted field into the metal can be seen. At late times a component of the field running along 
the metallic boundary is also visible and these may be related to surface plasmon polaritons [25]. 
No exact solution exists and this example is intended to show that the time stepping scheme can 
indeed handle a Drude model in a stable way. 
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Figure 4: Mesh for the Drude model (left) and Sellmeier model (right). 


4.2.2. Scattering from a Sellmeier dielectric 

The next model is an example for the Sellmeier model of a dielectric (see (23)) taken from [26] 
where a circular scatterer of constant permittivity sits on a glass substrate. The dimension of the 
unit cell and other parameters are chosen so that Fourier components of an incident wave with 
wavelength greater than 650 nm are scattered predominately in a different direction to those with 
a wavelength below 650 nm and the device is called a spectrum splitter. 

The geometric configuration is shown in the right panel of Fig. 4 where the parameters are as follows: 
L — 560 nm and R — 168 nm (so R/L = 0.3). In a slight departure from the optimized result in 
[26] we introduce a 20nm gap between the circular scatterer and the gas substrate (otherwise gmsh 
could not generate a mesh). The mesh used for the upcoming results is shown in the right hand 
panel of Fig. 4. 

We set n\ — 1.8, ri 2 = 1, c = 0.3 ^m/femtosecond and T — 8 femtoseconds and using 512 time 
steps so At = 0.0156. We chose 10 modes for the top and bottom boundaries. The glass substrate 
is assumed to be SF11 glass [27] so the Sellmeier model is 

2 , 1.73759695A 2 0.313747346A 2 1.89878101A 2 

6 7i~ = 1 —U - -1— - -I- - 

r 9 A 2 - 0.013188707 A 2 - 0.0623068142 A 2 - 155.23629 
with A denoting the free space wavelength in units of /^meters. More generally the Sellmeier model 
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j =1 

Converting first to the Fourier frequency domain and then to the Laplace transform domain this 
becomes 


6 r — 


1+ E 


Z - J 1 I Pj r>2 

j= 1 1 + 47r 2 C 2 S 

As explained earlier, the above rational function easily fits into our theory. 
We use the following function as the incident wave: 

f(r) = sin(2.899r)e" 2(T " 3)2 , r6l, 


(24) 


and the angle of incidence is 6 — 6°. For a graph of this function see Fig. 6. Considering the Fourier 
transform of this function, it has a maximum Fourier coefficient at the switch point 650 nm and 
hence has Fourier components on either side of the switch. 

Results are shown in Fig. 7. Again we do not have an exact solution, but demonstrate that the 
solution is stable. Spectral splitting [26] may be visible as two higher intensity zones moving in 
different directions (see t — 5.29,6.14 (a high intensity region moving up and to the right) and 
t — 6.99, 7.84 (a lower intensity region in the top left moving up slightly to the left)). 


5. Conclusion 

In this paper we have proved a general existence and continuous dependence result for time domain 
solutions of the grating problem with frequency dependent coefficients. We have shown that this 
leads to a convergent time discretization. Following spatial discretization with a typical non¬ 
overlapping finite element and spectral domain decomposition technique we have verified the time 
stepping convergence rate, as well as demonstrating the apparently stable numerical solution of two 
problems involving typical frequency dependent material coefficients. 

The study needs to be completed by an analysis of spatial discretization including mesh truncation 
and this will be the subject of a future paper. Further testing of the approach on real solar voltaic 
devices would also reveal the benefits and limitations of the time domain approach. In particular we 
need to investigate the relative efficiency of time domain and frequency domain approaches. 
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Figure 5: Snapshots of density plots the absolute value of the total field u;(x, t) for the Drude model 
at times t — 1.27,1.72,2.17,2.62,3.07,3.52 (fro® top left to bottom right). At later times faint 
waves traveling close to the surface of the metal are suggestive of surface plasmon polaritons [25]. 
































































Figure 6: A graph of the function /(r) against r for the incident field given in equation (24). The 
incident wave w l is given by replacing r by t — Ld\jc — dzy/c as shown in (4). 
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Figure 7: Snapshots of density plots of the absolute value of the total field w(x,t) for the Sellmeier 
model at times t = 3.59,4.44, 5.29, 6.14, 6.99, 7.84 (from top left to bottom right). 


























































































































